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Abstract 

Motivated by the well-known Papoulis-Gerchberg algorithm, an iter- 
ative thresholding algorithm for recovery of sparse signals from few ob- 
servations is proposed. The sequence of iterates turns out to be similar 
to that of the thresholded Landweber iterations, although not the same. 
The performance of the proposed algorithm is experimentally evaluated 
and compared to other state-of-the-art methods. 

1 Introduction and Problem Definition 

In many real world problems, instead of the complete signal, we have some 
observations of the signal of interest from which we want to reconstruct the 
original signal. In the simplest case, which is fortunately applicable to many 
practical situations, the observation process can be approximated as a linear 
operator: 

o = Kf (1) 

where / is the original signal of interest, o is the observed features, i.e. samples, 
and K is the (linear) observation operator. We are interested in the case where 
the number of observations is much fewer than the length of the original signal. 
That is, K mxn has much fewer rows than columns, i.e. m <C n. Furthermore, 
in practice the observation process is not exact and typically we can only obtain 



a distorted version of the observed features. This distortion is usually modeled 
by an additive error term. So the observation process can be modeled as 



g = o + e = Kf + e 



(2) 



in which e represents the (additive) observation error, e.g. noise. 

The objective is to find the original signal, /, from the set of available 
observations, g, while K is also known. That is, to solve the linear inverse 
problem ([2]). In order to do so one might minimize the discrepancy A(/) = 
|| Kf — g\\ 2 . However, except for the very special case that the operator K has 
a trivial null space (see |2 for details), the minimizer is not unique. In order to 
address this problem, i.e. to regularize the inverse problem, one might suppose 
a priori assumptions and impose different constraints on the solution, which is 
usually taken into account by adding a penalization term to the discrepancy. 
In this letter we are especially interested in the case where we have a sparsity 
constraint on the solution. For the sake of a more precise explanation, suppose 
there exists an orthonormal basis (?/>») in which the signal of interest can be 
expanded with the expansion coefficients =< f,ipi >; a large number of 
which are zero or negligibly small. In order to make use of this constraint for 
solving the inverse problem ([2]|, define the P-norm ||0|| p = (^,- |#i| p ) 1 ' p . One 
might then find the minimizer of the following functional, as the solution of Q 
with the aforementioned sparsity constraint: 



where < p < 2 and \i is the regularization parameter that can be chosen based 
on application. 

This is a well-known problem and different approaches to solving it have 



= \\Kf-g\\ 2 + v\\9\\ 



p 



(3) 
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been proposed. Here we will not go through the details of the problem, which 
have been widely studied by other authors. The reader is referred to [2] for a 
comprehensive discussion of the problem. For the sake of consistency, we follow 
the same mathematical notations as those used in [2J. Furthermore, it is worth 
mentioning that beside the Zp-norm introduced above, some authors have as well 
used the total variation (TV) norm as the constraint. See, for example, |13j . 

There are several methods of recovery available in the literature, among 
of which iterative thresholding algorithms are an important class. Iterative 
thresholding algorithms are, more or less, based on a thresholded version of the 
Landweber iterations (see, for example, jl]). I.e. the sequence of iterates has 
the general form 

f n = S 1 (f n - 1 +K*(g-Kf n - 1 )) (4) 

where K* denotes the conjugate of K, and 5 7 is a thresholding operator. In 
[5J the authors prove the convergence of the above iterative algorithm to the 
(unique) minimizer of ([3|, when S* 7 is the soft thresholding operator (see the 
definition of soft thresholding in section [2]) . Noticeable effort has been put into 
accelerating the original algorithm. In [5] the authors propose a method for 
accelerating thresholded Landweber iterations, which is based on alternating 
subspace corrections. Other methods for this purpose are introduced in [B], and 
[7 . Although use of soft thresholding is more common, some authors have, as 
well, used hard thresholding to address the above inverse problem. See [5], [S] 
or [TU] as examples. 

2 Description of the proposed algorithm 

In order to explain the underlying idea of the proposed method, let us begin 
with the following problem. In jl] Papoulis introduces an iteration method for 
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reconstruction of a band- limited signal from a known segment. Suppose /(f) 
is a signal of which we only know a small segment, g(t) = p T (t)f(t), where 

f 1 |*i < T 

Prif) — \ ■ Also suppose F{ijj) is the Fourier transform of /(f) 

otherwise 

and F(oj) = for \uj\ > a (bandlimitedness) . The objective is to reconstruct 
/(f) from g(t). 

In order to solve this problem, we begin with G(w), the Fourier transform 
of 17(f), and form Hi(uS) = G(cu)p a (uj), i.e. truncate G(oj) for \ui\ > a. In 
other words, we change g so that it satisfies the constraint on the original signal 
(bandlimitedness in this case). /11(f), the inverse transform of iJi(w), is then 
used to form /i(f) = g(t) + h\{i) — hip T (t), which recovers the known segment 
°f /■ /i(f) is supposed to be a better estimate of the desired signal, /(f), than 
g(t). This estimate can be further improved by repeating the above procedure 
in an iterative manner. That is, in the nth iteration, we form the function: 

H n (u) = F„_i(w)p (T (w) (5) 

compute its inverse transform, h n (t), and recover the known segment of the 
original signal 

/«(*) - g(t) + M*) - M*)Pr(f) (6) 

It can be proved that F n (ui) tends to F(ui) as n — ► 00 p]. 

In brief, in each iteration, we change the latest estimate of the desired sig- 
nal, i.e. the output of the previous iteration, so that it satisfies the constraint 
(bandlimitedness in this case). Since this process might affect the entire sig- 
nal, including the known segment, the known segment is then recovered before 
further progress. 

This problem is obviously different from our original problem stated in ([3]), 
because, firstly, it concentrates on the special case of recovering a continuous 



4 



signal from a known segment and, secondly, the constraint on the signal is 
bandlimitedness while the constraint of ^ is sparsity. Nevertheless, we will 
implement the above idea to solve our own problem as explained below. 

Based on the above algorithm, our iterative algorithm involves two main 
operations in each iteration, namely, an operation to maintain the constraint 
followed by an operation to recover the original observations. Since we are 
interested in problems with sparsity constraint, a thresholding operation can 
maintain this constraint for us, i.e. 



(7) 



where /" 1 is the latest estimate of the original signal, obtained in the previous 
iteration, and Sj is the soft thresholding operator, defined as 



(8) 



X + g X < 2 







W<3 



where s 7 (x) = < 

Analogous to (|6|), the original observations are then recovered by 



/" = h n + K*(g-Kh n ) 



(9) 



The sequence of iterates can, thus, be expressed in the following form: 



r = K*(g + S 1 (f n - i )-KS 1 (f n - 1 )) 



(10) 



with f° = K*g. 



Although (10 1 is not exactly a sequence of Landweber iterations, it can still 
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be viewed as a modified version of the thresholded Landweber iterations. Note, 
especially, the analogy between (10) and Q. 

In this letter we only introduce the algorithm and experimentally evaluate 
its performance, compared to similar state-of-the-art algorithms. A detailed 
discussion of the convergence of the iterative algorithm and its relation to the 
thresholded Landweber iterations is beyond the scope of the current letter and 
will be postponed to future publications. The motivation behind the proposed 
algorithm was briefly discussed, though. 



3 Experiments 

In all the experiments described below, the thresholding operator is applied 
to stationary wavelet transform (SWT) [11 coefficients, obtained using DB1 
(Haar) mother function for 1 level of decomposition. All thresholds are obtained 
using the well-known Birge-Massart strategy [12]. The iterative algorithm con- 
tinues until a convergence criterion, e.g. ||x fe+1 — sc ||/||a; || < i5, is met. For the 
sake of comparison, the results are compared with those obtained by L\ norm 
minimization [3] and total- variation (TV) norm minimization |13j . which are 
two well-known state-of-the-art methods of sparse signal recovery. 

Due to space constraints, the results of the experiments are included very 
concisely. More comprehensive results can be found at http : //mkayvan . googlepages . 
com/ sparsesignalrecovery 

3.1 Recovery of ID signals 

First, we consider the ideal case of sampling with no distortion, i.e. we assume 
e = in (J2j . The HeaviSine test signal (figure [T]) , from the well-known Donoho- 
Johnstone [2] collection of synthetic test signals, is reconstructed from different 
numbers, M, of randomly selected samples. Table[l]shows the the mean squared 
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Figure 1: Left: HeaviSine test signal; right: Shepp-Logan Phantom 



^Sample 


Rcc. by ( 


10 


) LI norm T.V. norm 


M = 70 
M = 100 
M = 150 
M = 200 


0.0339 0.13 0.0547 
0.055871 0.096 0.024 
6.51e- 03 0.038 0.0131 
4.65e-03 0.03 0.01 



Table 1: MSE between the reconstructed and the original signal for reconstruc- 



tion by (10) as well as by LI and TV norm minimization, from M observed 
samples. The original HeaviSine test signal constitutes N = 1024 samples. 



error (MSE) between the reconstructed and the original signal for reconstruction 
by (jlOj) as well as by LI and TV norm minimization. As it is obvious from 



results, reconstruction by (10 1 outperforms the two other methods in almost all 
cases. 



3.2 Recovery of 2D Signals 

Here we consider the classical problem of reconstruction of images from highly 
incomplete frequency domain observations, which has received considerable at- 
tention, because of its important applications in medical imaging, especially in 
MRI. In particular, acquiring MR images involves acquisition of 2-D Fourier 
domain data of the image. Due to the physics of the imaging device, this pro- 
cess is usually carried out by taking 1-dimensional slices from the 2-dimcnsional 
Fourier domain data of the image. This process is often too time-consuming, 
though, so for a rapid MR imaging it is desirable to take only a subset of these 
slices, e.g. a reduced number of Fourier domain samples taken over radial lines. 
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^Radial Lines 


Rec. by (10) T.V. norm LI norm 


K = 9 
K = 11 
# = 15 
K = 21 


24.9746 14.36 11.8 
29.2307 21 13.45 
39.3145 21.66 15.33 
199.7471 113.2541 27.1 



Table 2: Reconstruction of the 256 x 256 Shepp-Logan Phantom from samples 
along K radial lines in the 2D-DFT domain. Reconstruction quality is measured 
in terms of PSNR of the reconstructed image. 



Radial sampling is a common way of sampling over the 2-dimensional Fourier 
domain, and several authors have addressed this problem, especially as an im- 
portant application of compressive sampling. For the sake of comparison of the 
proposed method with other state-of-the-art reconstruction methods, here we 
will also address the same problem. The reconstruction problem might, how- 
ever, seem a bit different here, since sampling is confined to radial lines, instead 
of being random. Nevertheless, what we are essentially doing is reconstruction 
of the signal from its projections onto a lower dimensional subspace, which is 
exactly what was being done in the previous cases. 

The test image used is the Shepp-Logan phantom (figure [TJ of size 256 x 256. 
The pixels in this image take values between and 1, and the image has a 
nonzero gradient at 2184 pixels. The setup of our experiments is the same 
as that of [T3] which has been as well adopted by several other authors as a 
framework to evaluate the performance of their methods, including |15j . |16j . 
[T7], and [IB]. 

Table [2] shows the PSNR (peak signal-to-noise ratio) values of the recon- 
structed images from samples taken along k = 9, 11, 15, and 21 radial lines 
in the 2-dimensional discrete Fourier Transform (DFT) domain, compared to 
those obtained by LI and TV norm minimization. The reconstructed images 



by (10) as well as by minimizing 11 and TV norms, from samples taken over 9 



radial lines, are shown in Figure [2] 
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Figure 2: Reconstruction of the 256 x 256 Shepp-Logan phantom from sam- 
ples taken along K — 9 radial lines in the Fourier domain. From left to right: 



(a) Original 256 x 256 Shepp-Logan phantom (b) Reconstruction by (10) (c) 
Reconstruction by TV norm minimization (d) Reconstruction by 11 norm min- 
imization. 



dB 


Noisy image Rec. by ( 


10 


1 LI norm T.V. norm 


20 
30 
40 
50 


20.0663 26.2296 18.65 18.01 
30.0176 34.8228 22.94 25.89 
40.0028 43.7792 25.43 32.34 
50.0185 53.4432 26.15. 48.8 



Table 3: Reconstruction of the 256 x 256 Shepp-Logan Phantom from noisy 
samples affected by AWGN. The first column shows the PSNR of the AWGN; 
the second column shows the PSNR of the noise-affected image, from which the 
samples are taken, with respect to the original image. The corresponding values 
of PSNR, for reconstruction by (10 1, LI norm minimization, and TV norm 



minimization are shown in the third, fourth, and fifth columns, respectively. 

In the next experiment we consider the more realistic case of distorted obser- 
vations. Particularly, samples are taken from noisy versions of the Shepp-Logan 
Phantom test image, affected by additive white Gaussian noise (AWGN). The 
results are shown in table [3] Note that especially when the noise is significant, 
the reconstructed image enjoys considerably less error than the noisy one from 
which we got our samples. 



4 Conclusion 

Motivated by the Papoulis-Gerchberg algorithm, a method for recovery of sparse 
signals from very limited numbers of observations was proposed. Iterative 
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thresholding algorithms have been widely used to address this problem. Our al- 
gorithm also takes advantage of thresholding to maintain the sparsity constraint 
in each iteration. The signal is then reconstructed by iteratively going through 
a constraint-maintaining operation followed by recovery of the known features. 
The performance of the method was experimentally evaluated and compared to 
other state-of-the-art methods. 
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